
#define _THRESvb 0.0000001

class VB
{
	public:

	VB(void)
	{	
		boost::math::normal_distribution<> *d=new boost::math::normal_distribution<>(0,1); 
		_d=d;
		_conv=0;
	}
	~VB()
	{
		delete _d;
	}
	void operator()(mat X,double *Y,mat m,double si)
	{
		double e=1;
		double thres=_THRESvb;
		int p=X.n_cols;
		int n=X.n_rows;
		mat v(p,p);
		v.fill(si);
		v=diagmat(v);
		mat theta(p,1);
		mat Ym=vect2mat<double>(Y,X.n_rows,1);
		theta=inv(inv(v)+X.t()*X)*X.t()*Ym;
		mat b(p,1);
		int j=0;
		mat z;
	
		while(e>thres & MAXIT>j)
		{
			z=X*theta;
			for(int i=0;i<n;i++){
				if(Y[i]==1)
				{
					z(i,0)=z(i,0)+pdf(*_d,z(i,0))/(1-Phi(-z(i,0)));
				}else{
					z(i,0)=z(i,0)-pdf(*_d,z(i,0))/(Phi(-z(i,0)));

				}
			}
			b=inv(inv(v)+X.t()*X)*X.t()*z;
			e=abserr(theta-b);	
			theta=b;
			j++;
		}
		if(j>MAXIT){cout << "reached max iteration!";_conv=1;	}
		_m=theta;
		_S=-inv(Hessian(X,Y,theta,si));
	cout << _m;	

	}

	inline double abserr(mat X)
	{
		double sum=0;
		int n=X.n_rows;
		for(int i=0;i<n;i++)
		{
			sum+=abs<double>(X(i,0));
		}
		return sum;
	}

	inline double Phi(double x){

			double t=cdf((*_d),x);
			if(t==1)
			{
				t=0.9999999;
			}
			if(t==0)
			{
				t=0.00000001;
			}
			return t;
	}

	inline	mat Hessian(mat X,double *Y,mat theta,double si)
	{
		int p=X.n_cols;
		mat sum(p,p);
		sum.fill(0);
		int n=X.n_rows;
		double xb=0;
		for(int i=0;i<n;i++)
		{
			
			xb=as_scalar(X(i,span::all)*theta);
			double pxb=pdf(*_d,xb);
			double gpxb=cdf(*_d,xb);
			mat xx=X(i,span::all).t()*X(i,span::all);	
			if(Y[i]==1)
			{	
				mat temp=(xb*pxb/gpxb+pow(pxb,2)/pow(gpxb,2))*xx;
				sum=sum+temp;				
			}else{
				mat temp=(-xb*pxb/(1-gpxb)+pow(pxb,2)/pow(1-gpxb,2))*xx;
				sum=sum+temp;				
			}
		}
		return -sum-1/si;

	}
	mat Get_Sig(void){return _S;}
	mat Get_Mu(void){return _m;}
	int Get_conv(void){return _conv;}
	private:
	mat _m;
	mat _S;
	int _conv;
	boost::math::normal_distribution<> *_d; 
};

